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ABSTRACT 

Observations of disk galaxies at z ~ 2 have demonstrated that turbulence driven by gravitational 
instability can dominate the energetics of the disk. We present a ID simulation code, which we have 
made publicly available, that economically evolves these galaxies from z ~ 2 to z ~ on a single 
CPU in a matter of minutes, tracking column density, metallicity, and velocity dispersions of gaseous 
and multiple stellar components. We include an H2 regulated star formation law and the effects of 
stellar heating by transient spiral structure. We use this code to demonstrate a possible explanation 
for the existence of a thin and thick disk stellar population and the age- velocity dispersion correlation 
of stars in the solar neighborhood: the high velocity dispersion of gas in disks at z ~ 2 decreases along 
with the cosmological accretion rate, while at lower redshift, the dynamically colder gas forms the low 
velocity dispersion stars of the thin disk. 

Subject headings: galaxies: evolution - galaxies: ISM - instabilities - ISM: kinematics and dynamics - 
turbulence 



1. INTRODUCTION 

In the past decade, observations of galaxies near z ~ 2 
have revealed compelling evidence for the importance 
of gravitational instability in their dynamics and evo- 
lution. A whole class of galaxies has been observed 
who se images are dominated by large luminous clumps of 
gas (lElmegreen et al.ll2004l2005HForster Schreiber et all 
l2009f ). while measurements of the velocity dispersion of 
such massive star-forming galaxies have found values 
near 50 km/s spread acr oss the entire disk (|Cresci et alJ 
12001 iGenzel et al.l[20ll . This is difficult to reproduce 
with supernova feedback, which is strongest near the cen- 
ters of galaxies where the star formation rate peaks, and 
which is only strong enough to d rive velocity dispersions 
of ~ 10 km/s (poung et aLll2009l ) Other forms of stellar 
feedback may drive turbul ence ([Thompson et all 120051 : 
Elmc green fc Burkerdl201Ch . but we will concentrate on 
the case where turbulence is driven by gravitational in- 
stability in the disk. 

To a first approximation, the gravitational stability of 
a thin disk to axisymmetric perturbations is described 
by Toomre's Q parameter Q = Kcr/(7rGE), where k is 
the epicyclic frequency, a is the Id velocity dispersion, 
and S is the gas surface density. The disk is unstable 
when Q < 1. The importance of gravitational instability 
in high redshift galaxies arises from the high cosmolog- 
ical accreti on rates they exp erience, which drive up the 
value of E (jDekel et alj|2009f ). This instability gives rise 
to clumps of the sort observed at high redshift. The in- 
homogenous and time-varying gravitational field drives 
turbulence throughout the disk, regardless of the stellar 
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density or supernova rate. The energy source for these 
random motions must ultimately be the gravitational po- 
tential of the galaxy, so gas is transported inwards. 

Co smological simulations with sufficiently high resolu- 
tion (jBournaud fc Elmegreenll2009l : [Ceverino et al.ll2010f) 
successfully reproduce disks in which gravitational in- 
stability forms clumps and causes the inward migra- 
tion of material th rough galactic disks. Simulations 
of isol ated galaxies (jBournaud et al.l 1201 U iDobbs et al.l 
1201 lal lbT) with initial conditions set such that Q < 1 pro- 
vide a higher resolution view of such galaxies over a few 
outer orbits of the disk. These studies, while illuminat- 
ing, are expensive, since they must solve the equations 
of hydrodynamics in three dimensions over cosmologi- 
cal times. The model we present here solves the hydro- 
dynamical equations in the limit of a thin axisymmet- 
ric disk. Since quantities vary only in the radial direc- 
tion, the problem is computationally much cheaper to 
solve, allowing us to explore parameter space efficiently, 
while still solving the full ID equations of fluid dynam- 
ics instead of relying entirely on semi -analytic models 
(jDekel et al.ll2009t ICacciato et alj|201ll) 

Past ID models of gravitational instability in disks 
have a number of shortcomings. The rate at which mass 
and angular momentum are transported inwards is often 
parameterized and fit to the results of hydrodynamical 
simulations, rather than being derived from first princi- 
ples. The rotation curves are only allowed to be either 
Keplerian or flat. Energy is frequently assumed to be 
instantaneously equilibrated, which neglects the possi- 
bility that it might be advected through the disk. The 
pressure support of the disk is often treated as coming 
from thermal pressure rather than supersonic turbulence. 
Few models take into account the stellar component of 
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the disk, which becomes increasingly important as the 
galaxy evolves towards the present day, and can ulti- 
mately provide a variety of observable predictions. 

In particular, the age-velocity dispersion-mctallicity 
correlation of st ars in the solar neighborhood 
(jNordstrom et all 120041 ) . might well be explained 
by m eans of gravitational instability in high redshift 
disks (|Bournaud et al.l [20090 . The high velocity disper- 
sion in these disks means that the population of stars 
formed in that epoch wi l l star t with a high velocity 
dispersion (jBurkert et al.l H992). The gas disk cools 
as a result of slowing cosmological accretion rates, so 
younger stars are formed in a thinner, more metal-rich 
disk. This mechanism of thick and thin disk formation 
contrasts with the more common story that various 
secular processes and minor mergers heat thin d isk stars 
i nto a thick disk (e.g. [Binnev fc Tremainelll987[ ). 

iKrumholz fc Burkertl (|2010f ) (hereafter KB10) found 
an analytic steady-state solution to the full equations of 
fluid dynamics in the thin disk limit under the assump- 
tion that the disk self-regulates to maintain Q = 1. To 
make the problem tractable analytically, however, they 
required a handful of simplifying assumptions: they use 
an analytic approximation to Q, which becomes progres- 
sively worse at lower redshift as the ratio of gas to stellar 
velocity dispersion deviates from unity. They also as- 
sume that the velocity dispersion of stars and gas are 
equal, and the gas fraction at all points in the disk re- 
mains constant in radius and time. In this paper, we re- 
lax these assumptions and include treatments of stellar 
migration, metallicity, the non-zero thermal temperature 
of the gas, and evolution of individual stellar popula- 
tions. These improvements along with an efficient simu- 
lation code allow us to realistically evolve disks from high 
redshift to the present day at minimal computational ex- 
pense. 

In section 2 we derive the equations governing the 
evolution of the gas over time. Section 3 presents 
the derivation of the equations governing the stellar 
dynamics. In section 4, we derive the evolution of 
metallicity in the gas and stars. In section 5 we dis- 
cuss how these differential equations are solved nu- 
merically, and in section 6 we present the results for 
fiducial parameters chosen to be similar to the Milky 
Way. We conclude in section 7. The code we de- 
scribe, named GIDGET for Gravitational Instability- 
Dominated Galaxy Evolution Tool, is available at 
http : //www .ucolick. org/~jf orbes/gidget .html 



2. GAS EVOLUTION EQUATIONS 

2.1. Basic Equations 

We first give a brief overview of the derivation of the 
evolution equations for the gas column density E and 
velocity dispersion a. For more details see KB10. The 
equations of mass, momentum, and energy conservation 
for a viscous star-forming fluid in a gravitational field 
are 

dp 
dt 

p^ = -VP- P V^ + V-T, (2) 



-V • (>v) - (f R + p)p. 
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where p, v, e, and P are the gas density, velocity, 
specific internal energy, and pressure respectively. The 
star-formation rate per unit volume at an Eulcrian point 
is p„, with a mass loading factor p equal to the ratio of 
gas ejected in galactic-scale winds to the star formation 
rate. We will be employing the instantaneous recycling 
approximation (see section 0]), which approximates all 
stellar evolution as occurring immediately. Of the mass 
which forms stars, the gas will only lose the so-called rem- 
nant fraction, to stars, while the remaining (1 — /r) 
will be immediately recycled into the ISM. The gravi- 
tational potential is ip, T is the viscous stress tensor, 
$ = T lk (dvi/dxk) is the rate of viscous dissipation, and 
r and A are the rates of radiative energy gain and loss 
per unit volume. 

For a thin disk, we formally have p = T,S(z) and v z = 0. 
By expanding the fluid equations in cylindrical coor- 
dinates, integrating over z, assuming axisymmetry and 
v r <C «0, and dropping time derivatives of the poten- 
tial and the circular velocity, we can obtain evolution 
equations for the gas column density and gas velocity 
dispersion. The evolution of column density is given by 
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where (3 = dlnv^/ dlnr is the power law index of the 
rotation curve, T = J 2Tir 2 T r ^dz is the viscous torque, 
Ef F is the star formation rate per unit area, and 
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is the mass flux. The second equality follows from the an- 
gular momentum equation, which is in turn derived from 
the (j> component of the Navier-Stokes equation (equation 
see KB10). 

The derivation of the velocity dispersion evolution 
equation requires an equation of state, which we take to 
be P = pa 2 . The velocity dispersion has a thermal and a 
turbulent component. It is a reasonable approximation 
to treat both as contributing to the pressure so long as 
we average over scales much larger than the character- 
istic size of the turbulent eddies, which will be of order 
the disk scale height. 

Taking the dot product of the velocity with the Navier- 
Stokes equation and adding it to the internal energy 
equation yields an equation for the total energy, i.e. 
internal energy, kinetic energy, and gravitational po- 
tential energy. By decomposing the velocity as v 2 = 
v r + v 4> + ^°nt > the kinetic plus thermal energy may be 
rewritten 



(6) 



where the velocity dispersion is taken to be the quadra- 
ture sum of a thermal and non-thermal component, 
a 2 = a 2 +a 2 t . Neglecting the v 2 term as small compared 
to both a 2 and in a thin, rotation-dominated, Q ~ 1 
disk, employing radial force balance to set dip/dr = vi/r, 
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assuming a constant potential to set dv^/dt = 0, and in- 
tegrating over z yields the evolution equation 
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To fully specify the evolution of the gas, we need to 
set a rotation curve, a prescription for radiative energy 
gain and loss per unit area, and a procedure for find- 
ing the viscous torque. This will allow us to specify 
v<t>, ft-, Q = fTdz, C = J Adz, and T. The rotation 
curve is specified at run-time, and T is set by our treat- 
ment of gravitational instability (see section I2.2[) . We 
set Q = 0, which is equivalent to requiring that the en- 
ergy balance in the gas is completely determined by the 
effects of the viscous torque and radiative loss. We as- 
sume that the loss rate, meanwhile, is proportional to 
the kinetic energy density per disk scale height cross- 
ing time, in agreement with the decay rate of turbulence 
observed in full 3D MHD simulations of superso nic tur- 
bulence (|Mac Low et alj|1998t IStone et al.lll998D . 
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where 77 is a free parameter of order unity. If the decay 
time is exactly the crossing time, 77 = 1.5, since the ki- 
netic energy surface density is (3/2) Ecr 2 . In dropping the 
time derivative of at , we have assumed that the thermal 
velocity dispersion is unaffected by radiative dissipation, 
i.e. that the gas is isothermal. 

The scale height H is approximated by taking the so- 
lution to the equations of vertical equilibrium for a single 
component disk, H\ — a 2 /irGY,, and adopting it to mul- 
tiple components: 



H = 



ttG(£ + /X» 



(9) 



where / represents the relative importance of the stellar 
mass, or the stellar mass within a gas scale height. Tak- 
ing / = a I a* interpolates between two extreme cases: 
when a /a* -C 1, the scale height should approach the 
single-component value, i.e. / = 0. When a ~ cr», the 
two-component disk behaves (at least in terms of vertical 
density) like a single fluid with surface density £ + £*, 
i.e. / = 1. Note that the stellar scale height, which does 
not directly affect the dynamics of the disk, is just taken 
to be the single component solution, 
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which is reasonable for the small values of f g found within 
the star-forming regions of the disk. In reality the verti- 
cal structure of a self-gravitating disk in a dark matter 
halo is not this simple. However, excluding the effects 
of dark matter introduces an error of only 13%, even 
in the dark-matter dominated regions of the outer disk 



(|Naravan fc Joel[2002f) . Given the uncertainty in 77, this 
approximation is adequate. 

Substituting for the scale height and a 2 lt 
obtain a radiative loss rate of 



a 2 -a 2 , we 
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In this form, the radiative loss rate is the gas kinetic en- 
ergy per dynamical time multiplied by a factor to account 
for the effect of stars on the disk's thickness and a factor 
to zero out the radiative losses when there is no turbu- 
lence. As the gas velocity dispersion falls towards the 
constant thermal velocity dispersion, non-thermal mo- 
tions die away, the gas no longer dissipates its energy via 
shocks, and C — > 0. The gas temperature used to calcu- 
late at is a free parameter of the model, but fiducially we 
assume T g = 7000 A', appropriate for the warm neutral 
medium of the Milky Way. At high redshift when the 
gas is virtually all molecular, T <C 7000AT, but in that 
regime at /a <C 1 anyway, even if we use the higher-than- 
appropriate gas temperature. The choice of a t therefore 
has virtually no effect on the high-redshift evolution of 
the disk. 

The governing equations for the gas (equations 0] and 
[7]) are derived under the assumption that v z — 0. We 
therefore implicitly neglect the gravitational potential 
energy of the disk associated with its vertical extent, and 
the associated P dV work that the gas performs when it 
changes its scale height. Qualitatively, the effect of in- 
cluding these terms would be to provide the gas with 
another place to store energy which it gains when falling 
down the galaxy's potential well, aside from turbulent 
motion. Thus with these effects the gas velocity dis- 
persion would be slightly lower, and hence so would the 
dissipation rate, the gas column density, and the star for- 
mation rate. The dissipation rate and star formation rate 
are each already controlled by a free parameter which is 
uncertain at the factor of two level, so we are content to 
neglect these additional repositories of energy. 

2.2. Gravitational Instability 

The stability against gravitational collapse of a self- 
gravitating disk is given by a Toomre Q-like parameter. 
Several such fragmentation conditions exist in the lit- 
erature. We a dopt a modifie d version of the condition 
determined by iRafikovl (|2001[ ) , wherein the stability of a 
multi-component disk is considered with the stars treated 
realistically as a collisionless fluid. 
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where i indexes an arbitrary number of stellar popula- 
tions, 4>i is the ratio of the ith stellar population's ve- 
locity dispersion to the gas velocity dispersion, Iq(x) is 
a modified Bessel function of the first kind, and the Q 
parameter for each component is defined by 
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ttGE,- 
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The epicyclic frequency is k — \/2(j3 + l)f2, and q = 
ka / 'k is the dimensionlcss wavenumber, where k is the 
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dimensional wavcnumber of the perturbation. Values of 
q, or cquivalently k, for which Q(q) < 1 are unstable for 
an infinitely thin disk, and the q which minimizes Q(q) 
corresponds to the least stable mode. It follows that if 
Qua/ = mm (Q(<7)) < lj the disk is formally unstable, 
while if Qu,af > 1, the disk is stable. 

Computing the value of Q requires a minimization with 
respect to q. Since Q and its partial derivatives must be 
calculated frequently (see equation 1161 below') . it is com- 
putationally expedient to use an approximate formula 
which does not require such a minimization. KB10 used 
Q _1 « Q w $ = Qg 1 + Q* 1 ' as proposed bv lWang fc Silkl 
(|1994D , but t his approximation become s inaccurate when 
er g /er* < 0.5. iRomeo fc Wiegertl (|2011[ ) have proposed a 
more accurate approximation 
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Q rp if Q*Tn > QgTg , 

w<> if Q t < O T ■ ^ ' 

(15) 



This formula includes corrections for the fact that the 
disk is not razor-thin, and T g . A disk of finite thick- 
ness is more stable against gravitational collapse because 
its mass is spread out vertically, so larger values of the 
Tj increase the value of Q for a gi ven set of column 
densit ies and velocity dispersions. IRomeo fc Wiegertl 
(|2011h give approximations to these correction factors, 
Tj w 0.8 + 0.7a z j/a r j. For simplicity we have assumed 
an isotropic velocity dispersion, so = T g = 1.5. Qrw 
and its partial derivatives are straightforward to compute 
and accurate over a wide range of a g j a* and Q g /Q*. The 
stability parameter as determined by Qn a f should also 
be modified to include the effects of disk thickness, so 
our code can use either Q ps Qn a fT or Q ps Qrw- 

Disks where gravitational instability dominates the dy- 
namics are exp ected to be self-regulated near Q = 1 
(jBurkert et all 120101 ) . A disk with Q < 1 develops in- 
homogeneities in the gravitational field, which exert ran- 
dom forces on gas in the disk, driving turbulence. The 
ultimate source of this energy is the gravitational poten- 
tial of the galaxy, so mass must move inwards. If the disk 
had Q<1, more mass would gather into inhomogeneities, 
thereby increasing the driving of turbulence, which sta- 
bilizes the disk, driving Q upwards. Meanwhile if Q<T, 
mass transport through the disk slows even if the cosmo- 
logical accretion rate does not, which tends to add mass 
and destabilize the disk. We therefore take as a hypoth- 
esis that Q is a constant of order unity at all points in 
the disk at all times. Thus we can set 
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The evolution of the gas state variables E and a, derived 
in the previous section, depends on the viscous torque 
and its radial derivatives, so we can recast equation (jl( 
in the form 



dQ _jd 2 T 
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f o T-F = 0, 



(17) 



where the fi are coefficients which can be read off from 
the gas evolution equations, and F encompasses all terms 
which do not depend on the viscous torque, including all 
stellar processes, discussed in the following section, and 
the rate of radiative dissipation. In particular, 

a dQ I dQ 
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Usually F will be dominated by the first term, the 
radiative dissipation of energy, which tends to destabilize 
the disk by "cooling" the gas, making F > 0. In this case, 
one can interpret equation (fTTj) as requiring the torques 
to move gas such that it stabilizes the disk to counter the 
effects of this cooling. 

Equation (|1T[) is a second order ODE requiring two 
boundary conditions. At the outer edge of the disk, we 
specify the accretion rate of gas onto the disk, M ext ac- 
cording to a pre-calculated accretion history, typically a 
fit to average accretion histories from cosmological simu- 
lations (|Bouche et al.|[20lol ). The torque is related to M 
through equation ((SJ, so by rearranging that equation, 
evaluating quantities at the outer radius, and requiring 
a particular M ex t, we obtain the outer boundary condi- 
tion 

^) =-M ext v 4> {R)(l + f3(R)). (22) 



Here R is a fixed outer radius of the disk. This con- 
dition implicitly assumes that all gas is accreted at the 
outer edge of the disk, which is not an unreasonable ap- 
proximation as long as gas accretes mostly through cold 
streams 

At the inner boundary, we require that the disk and 
bulge exert no torques on each other, 



ea =rn = o 



(23) 



The inner edge of the computational domain is ro, cho- 
sen for numerical reasons to be non-zero. Note that 
this boundary condition is somewhat different than the 

one used in KB10, namely (T) r =r = — M ex tV<f,{R){l + 
(3(R))tq for a flat rotation curve. This will approach the 
physically motivated value of equation (|23[) in the limit 
that ro — > 0, and was chosen to satisfy a regularity condi- 
tion at the inner boundary. However, since our goal here 
is not to obtain an analytic solution, there is no need 
to impose such a condition. In practice we have experi- 
mented with both choices in our numerical calculations, 
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and we find that the choice of inner boundary condition 
has negligible effects at radii r tq, which is the great 
majority of the disk. 



3. STELLAR EVOLUTION EQUATIONS 

In addition to the gas, we would like to know how 
stellar populations in the disk evolve with time. The 
stars will provide most of the observable consequences of 
the model, in addition to determining, along with the gas, 
whether the disk is gravitationally unstable. Among the 
questions we are interested in investigating is the cause 
of the age-velocity dispersion correlation, namely that 
older stars have higher velocity dispersions. Therefore it 
is useful to not only keep track of the stars as a single 
population with a single column density £* and velocity 
dispersion cr* (each a function of radius and time), but 
also to bin the stars by age, so that and cr*^ describe 
the ith age bin. 

The overall stellar population, along with each sub- 
population, will be directly affected by two processes - 
star formation and stellar migration. The two effects 
may be added together, recalling that of the gas which 
forms stars, only a fraction fn will remain in stars after 
stellar evolution has taken its course, 



alteration: 



SF 



^Mig 



(24) 



Evolution equations for each process will be derived sep- 
arately below. 



3.1. Star Formation 

The rate of star formation will depend on the proper- 
ties of the gas from which stars form. In particular, in 
a sufficiently large region of the disk, the star formation 
rate will be proportional to the molecular gas mass di- 
vided by the free fall time, defined to be y/3n/(32Gp). 
In deriving the gas evolution equations, we assumed that 
formally the density was given by Y>5(z), but this is of 
course an approximation. The disk will have a finite 
thickness of order the scale height (defined by equation 
so we take the density to be p = Y,/H. Thus we can 
write the star formation rate density 



Sf F = e ff / ff2 lV32Gp/(3^) = £ff foE K \/32/3 ( 1 + 

(25) 

For molecular gas, the effici ency of star formati o n per 
free-fall time is e ff ~0 . 01 ((Krumholz fc McKeel 120051: 
iKrumholz fe Tanl l2007t iKrumholz et al.1 120111 ). though 
this may be significantly higher or lower given ob- 
ser yational uncerta i nties. Foll o wing progress made 
bv IKrumholz et all (|2008l l2009al ). IMcKee fc Krumholzl 
have analytically approximated the molecular 
fraction of the gas, Jh 2 &s a function of metallicity and 
surface density. We adopt this prescription with a slight 
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if s < 388/203 



otherwise 



(26) 



ln(l + 0.6x + 0.01x 2 ) 



0.6t> 



X = 3 



1 1 + 3.1(Z/Z ) 



0.365 



4.1 



t c = 320 c (Z/Z Q )(E/1 g cm" 



(27) 

(28) 
(29) 

where Z is the gas metallicity. We take the solar metal- 
licity to be Zq — 0.02, and c encapsulates the effects 
of clumping in the gas when averaging over large re- 
gions. Since the model presented in this paper takes 
averages over large areas of the disk, w e take c ~ 5, as 
determine d in IKrumholz ~eT ah (l2009bft . The modifica- 
tion from IMcKee fc Krumholzl (t2O10F is that we impose 
a lower limit on fjj 2 of 3%, motivated by the observa- 
tion that even extremely low total gas surface density 
regions form stars at a rate consiste nt with a constant 
H 2 depletion time (jBigiel et al.H201lD . 

Equation (f2"5")) is used to update the stellar column den- 
sity, and it also enters into the gas column density equa- 
tion (equation [4j through the conservation of mass. At 
any particular time in a simulation, all but one of the 



0. Formally we can write this as 



^SF 



G(A(t) - Ay 



i) &(A i dtl - A(t)) (30) 



where Q(x) is a step function, one for x > and zero for 
x < 0, A(t) is the age that a star will be at redshift zero 
if it forms at time t after the beginning of the simulation, 
and A youn g t i and A a id.i are the boundaries of the ith age 
bin. 

To update the stellar velocity dispersion of a stellar 
population, we require that the new kinetic energy of 
the population be equal to the old kinetic energy plus 
the energy of the newly formed stars, 



_2 \ 

iJnew 



(31) 



where we have assumed that the newly formed stars have 
the same velocity dispersion as the gas from which they 
form. Setting £* ine ™ = S*,oid + fR{dT,f F ), we can rear- 
range, solve for er* ineu; , and expand to first order in the 
small quantity dT,^ F /S*. D ;, 
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'(SM<i)°M + /fl(<*S?,i> 

^*,i,old + f R (dZ s J) 
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/fl(rfSff) 



2S* i oldO * A old 



Thus in the limit that the time step and therefore the 
density of new stars produced is small, we may use the 
definition of a derivative to write 



SF 



h 



i 



2 XI sjt (7 sjt 



V-<7*i)£?,i for£.,i>° 

(33) 

We only need this derivative for its contribution to the 
torque equation (equation ITB)) . in which it will always be 
multiplied by the term dQjdu^^. To actually update the 
quantity cr*,,, we use the exact relation of equation (|3Tj) , 
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which holds even if E*^ = 0. Note that when E*^ = 0, 
this new population of stars will have no effect on the 
torque equation, since dQ/da*^ = 0, i.e. non-existent 
stars do not affect the stability of the disk. Thus equation 
([3"3"[) need only be employed when E* ^ > 0. 



3.2. Radial Migration 
star formation, stars are subject to 



In addition to 
radial migration. In particular, when < 2, tran- 
sient spiral arms form which attempt to stabilize the 
disk 
[19851 



(ISellwood fc Carlberd [l98i ICarlberg fc Sellwoodl 
Sellwood fc Binnevl 120021) 7 N-body simulations 



(jSellwood fc Carlberdll984l ) suggest that this heating is 
such that 



dQ 
dt 



Mig 



max 



Qi 



(34) 



where T m i g is the time scale in local orbital times over 
which this heating occurs, typically a few orbits, and 
Qiim is the value of Q* above which the stars are sta- 
ble to spiral perturbations. Equation ([3"4"[) assumes that 
this mechanism acts independently of the torques which 
act on the g clS 3j result of the axisymmetric insta- 
bility described in section 12.21 In z ~ 2 galaxies with 
morphologies dominated by clumps containing both gas 
and stars, one might expect the axisymmetric instabil- 
ity to affe ct both components e qually, as assumed in the 
models of iCacciato et al.l (|2011l) . However, it remains an 
open question whether these clumps are disrupted on a 
dynamical timescale by a stellar feedback process, just 
like giant molecular cloud s , their low-redshift analogues 
([Krumholz h Dekell l20Tot iGenel et al.1 12012I ). Even if 
clumps are long-lived, they cont ain a relatively smal l part 
of the total stellar population ((Murray et al.|[20To() . and 
thus their impact on stellar migration might be small. 
Moreover, in most realistic situations, the scale height of 
stars will be significantly greater than that of the gas, so 
an instability dominated by the gas will have little trac- 
tion on the stars. As long as a* jo is appreciably greater 
than unity, which it is in our fiducial model (section 
we expect this treatment to be reasonable. 

The time derivative of Q* may be re-expressed in terms 
of the time derivatives of E* and cr* using the definition 
of Q*, 



dQ, 
dt 



Mig 




-Q J_^ Ml9 

\ cr* dt 



dt 

1 <9£* Mi9 \ 



E* dt 



(35) 



The partial time derivatives on the right hand side will 
depend on the mean velocity of stars in the radial direc- 
tion, v*^, and so the forcing imposed by equation ([3~4"[) 
will yield an ordinary differential equation for v*^ r . The 
value of f IS then used to evolve E* and <r». 

This approach assumes a single bulk velocity of stars 
in the radial direction at ea ch radius, y*. r (r ). It has 
been well-demonstr a ted (e .g. iBird et ahl 120121 ) that the 
ISellwood fc Binnevl (|2002f) mechanism scatters stars in 
both directions, i.e. a star born at some galactocentric 
radius may end up with a guiding center radius multiple 
kpc away. There are additional scattering mechanisms, 



such as two-body scattering and the resonant overlap 
between spirals and the bar (|Minchev &; Famaevl 120101 : 
iBrunetti et al.ll2011[ ). which will also redistribute stellar 
angular momenta. Modeling this redistribution is critical 
in explaining the detailed properties of Milky Way stel- 
lar populations. However, there are no straightforward 
prescriptions to model all of these effects. We therefore 
ignore for now the effects of radial mixing and merely 
require conservation of mass and energy, and that the 
stars will stabilize themselves if they are subject to spi- 
ral instabilities. 

The evolution of E* and function of f IS de- 

termined by the continuity equations for mass and energy 
of the ith stellar population 



<9E, 



Mig 



i)t 



dV E* iV* r 
r 



= 



(36) 



d_ 

oi- 
id_ 

r dr 



Mig 



[rE^tw (y$ + Sal , + 2^)] = (37) 



Expanding the energy equation using the product rule 
and employing the mass equation to cancel terms leaves 

[(^ + 3<, + 2V0] + 

E*,^*,,.^ [(vl + 3al, t + 2V>)] =0 (38) 

The time derivatives of and ij) are zero by assumption, 
so expanding the surviving derivatives, setting d?p/dr = 
v^/r and dv^/dr = ftv^/r, and rearranging yields 



da*,i MiB 
dt 



'(1 + j8)«2 aci., 



3?'C7» 



Or 



(39) 



The corresponding equation for stellar column density 
follows immediately from the continuity equation: 



nr..,, .. or 



dt 



— E* 



Or 



' — v*, r — — E H , i iW*, r /r (40) 



di 



Substituting the transport equations into equation ([35 
and imposing equation ([34)1 yields 



2-Kr- 



v l (1 



3r 



1 <9er» 
cr* dr 



2nr dv* 



1 <9E* 
+ E* dr 

max(Q Hm 



v<f> 



Or 



TMigQ 




This is a first order ordinary differential equation (since 
at any particular time we treat all variables as functions 
of radius only), requiring a single boundary condition 
which we take to be w*, r (r = R) = 0, which means that 
no stars are allowed to migrate between the outer edge of 
the disk and the IGM. This boundary condition guaran- 
tees that the bulk velocity of stars in the radial direction 
will be inwards at all radii, which means this method 
does not conserve angular momentum; to compensate 
for a large mass of stars moving inwards, a small mass of 
stars would need to move outwards. The error we make 
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in conservation of total angular momentum is about 2% 
in the fiducial case. 

4. METALLICITY EVOLUTION 

4.1. Advection of Metals in Gas 

To describe the evolution of the metal content, we be- 
gin by defining E^, the surface density of metals, so lo- 
cally the mctallicity of the gas is Z = Y^zfE- The conti- 
nuity equation for T,z is 



d 



1 d 



= 7. — ttMz - E z 
at 2itr Or 



Sz 



(42) 



where E§ F = Tif F Z is the rate at which metals are in- 
corporated into newly formed stars, and Sz is a source 
term for metals injected into the gas by supernovae and 
AGB stars. Note that, in writing this equation, we ne- 
glect transport of metals through the disk by either tur- 
bulent diffusion or galactic fountains. The inward flux of 
metallic mass is 

Z &T 

M Z = MZ = 7——^--^-, (43) 

*W + P) or 

which follows from equation [5j The left hand side of 
equation [42] can be reexpressed in terms of Z by noting 
dY,z/dt = ZdE/dt + mz/dt. Equation ED then be- 



27rr(l + p) z v^ \ or Or Or 



dlar 
dr dr 



"(1 + 0) 



d 2 r 

dr 2 



- Ef F + Sz (44) 



Comparing this with the previously derived gas surface 
density evolution equation, we can cancel most of the 
terms on the right hand side with ZdT^/dt, leaving only 



dZ_ 
~dt 



1 



dlnZdT 



(P + l)rE^ dr dr 



E ' 



(45) 



Inflowing gas has some metallicity Zigm, which we fix 
at Ziqm = O-l^fTj for the entire simula tion. Simulations 
flShen et al.l |2011[ ) and observations (jAdelberger et al.l 
120051) suggest that the circum-galactic medium is en- 
riched to this degree as early as z = 3. 

4.2. Metal Production 

For simplicity, we adopt t he insta n taneo us recycling 
approximation, proposed by iTinslevI ([1980D . to specify 
Sz, the production rate of metals. First we recognize 
that metals are produced in supernovae and AGB stars. 
To a first approximation, we can assume that the life- 
times of stars that dominate metal production are much 
smaller than the timescales with which we are concerned 
in this paper, so metals enter the ISM at a rate propor- 
tional to the star formation rate. Not all gas which forms 
stars is returned to the ISM, since low-mass stars do not 
leave the main sequence in a Hubble time and even high- 
mass stars form remnants. Defining the remnant fraction 
fn as the fraction of gas forming stars which will end up 
not being returned to the ISM, the surface density of 



recycled gas appearing in the ISM is (1 — /^)Ef F . Su- 
pernovae and normal stellar evolution will enrich a small 
fraction of this gas, namely um, the yield. The surface 
density of metal production is therefore 



Sz = 2M/C(1 - 



^SF 



(46) 

Assuming a lChabrierl ()2005|) initial mass function and a 
coarse approximati on for the ultimate f ates o f stars as 
a function of mass, iKrumholz &: Dekel (|2011| ) compute 
fu = 0.46. Assuming in addition a production function, 
the fraction of a star's i nitial mass converted to a given 
element, from iMaederl (|1992D . they compute a yield of 
ijm = 0.054 for Solar metallicity stars. The effective 
yield may be somewhat smaller than this, since galactic 
winds driven by supernovae tend to eject material which 
is richer than average in metals. The factor of £ < 1 rep- 
resents the ratio of metallicity in the ISM to mctallicity 
of ejected material. We adopt ( = 1, corresponding to 
an assumption that the ejecta are well-mixed with the 
ISM. This value will i n principle depend on the mass of 
the galaxy considered (|Mac Low fc Ferraralll999[ ). How- 
ever, owing to the high resolution required, to date no 
simulation has reliably calculated the degree to which 
metal-rich gas is preferentially ejected. Changes in the 
exact value of £ roughly translate into the normalization 
of the metallicity distribution in the gas, so our fiducial 
value of £ was chosen to give reasonable values for this 
normalization. 

4.3. Diffusion of Metals 

The metallicity gradients produced when accounting 
only for metal production by stars and advection by in- 
flowing gas are far steeper than the observed gradient 
in the Milky Way. Metals are formed in proportion to 
the star formation rate, which tends to be high towards 
the center of the simulated galaxies. Meanwhile the in- 
flow of gas throughout the disk concentrates the metals 
even further. To explain the relatively small observed 
metallicity gradients, one must allow metals formed at 
small galactic radii to reach large radii. This may oc- 
cur either in the plane of the disk (diffusion) or out of 
the plane (galactic fountains) . By assuming a fixed value 
of Zigm — 0.1Z©, we have already implicitly assumed 
some sort of transport of metals from the galaxy into 
its surrounding medium. However, rather than model- 
ing this transport in any detail, let us consider only the 
diffusion of metals through the disk. 

In general, a diffusion equation will have the form 



d d 2 
—M z = D—M z 
Ot Or z 



(47) 



where D is the diffusion coefficient and Mz = 27rrArEZ 
is the gas-phase metal mass in a given cell. At an order 
of magnitude level D may be estimated by taking the 
typical velocity of gas in the disk, cr, and multiplying 
by the typical length scale of perturbations, namely the 
2d Jeans scale, a 2 /GT,. For simplicity we simply adopt 
D = kzv c /,(R)R where kz and D will be constant at 
every time and location in the disk. Numerically we take 
kz = 10~ 3 which is of the correct order of magnitude 
and yields a metallicity gradient of order 0.1 dex/kpc 
(see figure [2])., which is com parable to observed v alues in 
isolated spiral galaxies fe.g. iZaritskv et aL| [l994') 
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4.4. Metals Locked in Stars 

The metallicity of a given stellar population can be up- 
dated when new stars are added to it by again assuming 
instantaneous missing. The new metallicity is just an av- 
erage of the old metallicity and the metallicity of the gas, 
weighted by the surface density of the extant stellar pop- 
ulation and the newly formed population respectively. 



Z* 



old* 



f R Z(dZ 



*,2 / 



J *,i.new 



(48) 



Here, as in equation (|3T|) - dEff = Ef-fdi, is the surface 
density of stars formed in a given time step in a given 
stellar age bin i. 

Besides the formation of new stars, a given stellar pop- 
ulation is subject to migration through the disk, as dis- 
cussed in section 15721 Since the stars migrate through the 
disk with a mean velocity set by equation (|4ip , the metal- 
licity profile of a given population of stars evolves under 
a continuity equation for the metal mass, 



-rs z 



\Mig 



df^E* %Z* iV* r j 
r 







(49) 



Subtracting the continuity equation for total stellar mass 
(equation I36p. we obtain 



dZ*, 
dt 



Mig 



dZ* 



Or 



(50) 



for the evolution of stellar metallicity. Equations (|48p 
and ([50]) fully describe the evolution of the metallicity of 
the ith stellar population. Note that these equations ne- 
glect radial diffusion of stars, only ta king into account the 
mean velocity v*. r - Ra dial mixing (|Sellwood fc Binnevl 
l2002tlRoskar et aLll2011f ) is required to explain the spread 
of metallicities in stars at a fixed age and radius, and un- 
doubtedly leads to a shallower stellar metallicity gradient 
than what we obtain. 

5. NUMERICAL METHOD 

5.1. Computational Domain 

In deriving the gas evolution equations, we assumed 
the disk to be thin and axisymmetric. Thus the disk is 
described by variables which change only in radius and 
time. We therefore define a mesh of radial positions 
with a fixed number of points, n x , logarithmically spaced 
between the outer edge of the disk at a fixed radius R 
and a fixed inner cutoff r m i n , usually chosen to be r m i n = 
OMR. Explicitly, 



= R 



R 



l-(i-l)/(n x -l) 



(51) 



The highest spatial resolution is therefore given to the 
region with the shortest dynamical times. 

Time, tracked in units of the orbital period at radius R, 
begins at zero when the simulations are started, typically 
at z = 2, and reach a few tens of orbits at z = 0, depend- 
ing on the assumed radius and circular velocity. The size 
of the time steps are calculated by first determining all 
timescales defined by dividing each state variable at each 
position by its time derivative, picking out the minimum 
timcscalc, and multiplying it by a small number TOL. 
usually taken to be 10~ 4 . Larger values of TOL lead 



to numerical instabilities near the inner boundary, which 
is especially susceptible to such issues because the local 
dynamical timescale in the disk is fi _1 oc r for a flat 
rotation curve. 



At- 



■■ TOL x mini 

An), 



dZ/dt 
a 



(n), 



(r<), 



da/dt 
0.01 



in), 



<9£*/cV da*/dV in TOL 



(52) 



A maximum time step of 0.01 outer orbits is imposed 
to prevent systems extremely close to equilibrium from 
advancing too quickly. 

5.2. PDEs 

At each time step, the code solves the equations in 
non-dimcnsionalizcd form (sec appendix) in the following 
order. First, we solve equation (|4ip to determine u*. r at 
all radii. The equation is of the form % = /io«*,r + 
dv^^ r /dr with 



H 



h - 



max(Q| im - Q*, 0)-^ 

2-KrTMigQ* 

/3) _ 1 da* 
dr 



v l (1 



3r 



1 dS. 
£* dr 



(53) 
(54) 



The boundary condition specifics v*, r at the outer edge 
of the disk. Thus rewriting the radial derivative as a 
finite difference and employing a backwards Euler step, 
we can write an explicit update equation, 

v*An-i) « (55) 



1 - (n 



)ho(f 



which we solve itcratively by starting with the specified 
value of v* >r (r n x) = and moving inwards. 

Using the value of v*^ r along with the current values 
of the state variables, we calculate the coefficients of the 
torque equation (equation 1 1 T[) . To solve the resultant 
linear PDE, we employ a similar finite difference method, 
which approximates 



9% ^ %+i — %-i 
dr r t+ i - r t -i 

dr 2 



r i+l/2 



r i-l/2 



T, 



i+l 



H % — 17-1 



t+i 



(56) 



(57) 



Since we are using a logarithmically spaced grid, r i+1 / 2 — 
y/fifi+i- By plugging these approximations into the 
torque equation, the problem reduces to the inversion 
of a tridiagonal matrix. 

The forcing term in the torque equation, (|21[) generally 
acts to destabilize the disk, since its largest term comes 
from radiative cooling of the gas and cooler gas is more 
prone to gravitational collapse. The torque equation re- 
quires that the gravitational torques exactly counteract 
this effect to maintain dQ/dt = 0. However, in the event 
that the forcing term in the torque equation becomes 
negative as a result of stellar migration and a reduced 
rate of cosmological infall leading to C — s- 0, we set it to 
zero so that the gas is not forced to destabilize the disk. 
This in turn allows positive values of dQ/dt. We do not 
allow the forcing term to return to the value given by (|2ip 
until that value is again positive and Q has been allowed 
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to rise and then fall back down to Q = Qf . This allows 
the simulation to follow disks which stabilize at least tem- 
porarily, for example because of a lull in the cosmological 
accretion rate, and then return to a marginally unstable 
state. For the smoothed average cosmological accretion 
history used in our fiducial run, parts of the disk which 
stabilize remain that way because the accretion rate is 
monotonically decreasing. 

With T, dT/dr, d 2 T/dr 2 , and u» ir , we can now eval- 
uate the derivatives of the state variables. Where radial 
derivatives of the state variables or other quantities ap- 
pear in the evolution equations or the coefficients of the 
above differential equations, a minmod slope limiter is 
used to evaluate them. In particular, if L = (-A(fj) — 
A{r t - l ))/{r, l -r l - 1 )&nAR={A{r l+1 )-A{r l ))/{r l+l -r l ) 



dA 
dr 



< L if \L\ < \R\ and LR > 
(n) = { R if \L\ > \R\ and LR > 
I otherwise 



(58) 



where A is a stand-in for any quantity. This strongly 
suppresses noise on the scale of the mesh separation by 
zeroing out rapid variations in the derivatives. 

With the time derivatives calculated at each point, we 
simply take a forward Euler step to update the state 
variables, namely S, a, Z, £*, cr*, and for each age- 
binned stellar population, cr*,;, and Z*j. Typical 
runs have time steps limited by the rate of change of 
the gas state variables near the inner boundary of the 
disk where the dynamical timescalc is shortest. On a 
single processor, runs take about one day to complete 
if we nume rically evaluate Q(q) and its derivatives us- 
ing the full IRafikovl (|2001D formalism. We can shorten 
this by an order of ma gnitude by using the appr oxima- 
tion to Q suggested by iRomeo fc Wiegertl (|201lD . This 
approximation is much more efficient because Qrw and 
its partial derivatives may be calculated as functions of 
the state variables alone, without the need to minimize 
over a wavenumbcr or compute the partial derivatives 
dQ/d{S, cr, cr*,i} numerically as required by the full 
Rafikov Q. 

5.3. Initial Conditions 

By assuming a fiat rotation curve, fixed gas fraction, 
equal stellar and gas velocity dispersions, a simple ana- 
lytic approximation to Q, and ignoring stellar processes 
(formation and migration), KB10 were able to compute 
an equilibrium solution to the evolution equations. In 
particular, 



= J_ ( GM extfi 



1/3 



irGr 



P a GM e 



1/3 



'/ 



(59) 
(60) 



Here M e xt,o is the accretion rate of gas onto the outer 
edge of the disk at the start of the simulation, and f g is 
the gas fraction, assumed to be constant in radius. By 
assumption, cr* = a and £* = S(l — f g )/f g - 

If we relax the assumptions that the velocity disper- 
sions of both components are identical and Q = 1, add 
a factor to correct for finite disk thickness, but retain 



the approximate form of Q for an infinitely thin disk, 
Q^ 1 py Qg 1 + Q" 1 , we obtain a modified version of the 
equilibrium column density, 



T Vj, 



<Pofg 



GM, 



1/3 



ex 1.0 



Qf ^Gr fg 



1+1 



Vfg 



(61) 



where <fo> = cr^/er is a free parameter , T ss 1.5 is the 
thickness correction, and Q f is the fixed value to which 
Q will be set everywhere in the disk. To initialize the 
simulations, we use equations (|59l) and (jFTj) . We then 
adjust cr* = (poa keeping 0o fixed until Q = Qf exactly at 
each cell of the grid. Finally, we run the simulation with 
stellar processes turned off, i.e. e// = Qum = 0, and with 
M ext fixed to its initial value, M ext fi, to allow the gas 
to adjust to an equilibrium configuration. The greatest 
effect of this adjustment occurs at the inner edge of the 
disk, since these relations were derived using a different 
inner boundary condition and under a more stringent set 
of assumptions. Once the state variables are changing 
sufficiently slowly, we have found our initial conditions 
and therefore return e//, Qum, and M ext (t) to their user- 
specified values. 

6. FIDUCIAL MODEL 

While our code is quite general, here we describe a 
simple model run using it in order to demonstrate its ca- 
pabilities. In future work we will explore a much wider 
part of parameter space, using more realistic cosmologi- 
cal accretion histories. 

6.1. Setup 

The formalism presented here requires a rotation curve, 
accretion history, and fixed inner and outer radii to be 
specified before the simulation is run. Since we employ 
a logarithmic computational grid, there is little cost to 
extending the outer radius out to 20 (as opposed to 10) 
kpc. This allows us to follow the transition of the outer 
disk from somewhat molecular at high redshift to atomic 
at low redshift. For the inner truncation radius, we take 
ro = O.Oli? = 200pc. The exact value will affect the 
quantitative results within a few kpc of the center of the 
disk, but the exact results of the simulation in this re- 
gion should be taken with a grain of salt anyway. Here 
cr* reaches a similar order of magnitude as the circu- 
lar velocity, which we take to be independent of radius, 
v <t>{ r ) = 220 km/s, so our treatment of this region as 
a thin disk is not valid. Moreover, the inner boundary 
value for the torque equation, which we take to be zero 
- no torque is exerted by the region within the trunca- 
tion radius on the disk - could easily be some small but 
non-zero value. 

The accreti o n hist ory employs the fitting formula from 
iBouche et ail (|2010f ). namely 



M(t) = 7 e m / b , .i8 Mfo (1 + z) 4 -" M Q /yr (62) 

where Mh, 12 is the halo mass in 1O 12 M0, fb,o.is is the 
baryon fraction of the accreting matter normalized to 
18%, and ti n is zero for Mh,i2 > 1-5 but varies linearly 
in time from 0.7 down to 0.35 between redshift 2.2 and 
1. Before redshift 2.2, ej„ = 0.7, and after redshift 1, 
ei n = 0.35. We choose /&,o.is = 1> and an initial halo 
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mass which will grow to be about 10 12 Mq at redshift 
zero. The formula governing the growth of the halo mass 
is given in the same paper, 

M h = 34.0 Ml'll (! + z f A M °/V r i ( 63 ) 

so an initial halo mass of Mu.n = 0.27 at z = 2 pro- 
duces a Milky Way-analogue galaxy with M/> 12 « 1 at 
z = 0. We note that some of the baryonic accretion may 
go into expanding the outer radius of the disk, instead 
of being transported inward, which would reduce the ac- 
cretion rate below the estimate given in equation (|62|) . 
However, since the baryonic mass of galactic disks out- 
side 20 kpc is generally a negligible fraction of the total, 
this clearly cannot be a large effect, and the error we 
make by neglecting it is small compared to the general 
uncertainty in the cosmological accretion rate. 

In addition to these functions, there are several free 
parameters controlling various physical processes in the 
disk. The star formation efficiency per freefall time is 
Efj = 0.01. The mass loading factor of winds ejected 
from the galaxy in proportion to the star formation rate 
is ij,= 1, c hosen to roughly correspond to observations 
(jErbl I2008D . The fraction of turbulent energy in the 
gas which will decay in a scale height crossing time is 
77/I.5 = 1. The time scale for a Q* = Qu m — 1 popu- 
lation to approach Q* = Qu m is T rnig = 2 local orbital 
periods, and the value of Q* below which the stars are 
subject to transient spiral instabilities is Qu m = 2.5. For 
computational convenience, we use Q w Qrw to evalu- 
ate the disk's stability. We will explore the sensitivity 
of the results to these parameter choices in future work. 
Here our goal is merely to demonstrate the method and 
its results. 

The value of Q everywhere in the disk is fixed to Qf. 
Theoretically Q is expected to be self-regulated to a value 
of order unity. Formal stability criteria derived from the 
perturbed equations of motion for infinitely thin disks 
find the disks to be unstable when Q < 1, so the marginal 
stability which we ass ume here would im ply Q = 1 . How- 
ever, recent work by lElmegreenl (|2011f) has shown that 
for a realistically thick disk where the gas cools on the 
order of a dynamical time, a marginally stable value of 
Q is closer to 2 or 3. This i s consistent with the obse rva- 
tional evidence compiled by [Romeo fc Wieeertl (I20T1 for 
nearby spiral galaxies, namely that when Qrw for these 
disks is calculated, the values typically fall between 2 and 
3 for most galaxies at most radii. Thus we adopt Qf = 2 
as a fiducial value. 

Finally, to specify the initial conditions fully, one must 
choose an initial gas fraction and a ratio of stellar to 
gas velocity dispersion. Since the only way the stellar 
velocity dispersion can decrease is by mixing it with a 
lower-velocity dispersion population, it is reasonable to 
expect this ratio to be greater than unity. The simplified 
models of gravitati onally unstable galaxi es evolving from 
z >• 1 discussed in lCacciato et al.1 ()2011[ ) suggest that by 
z ~ 2, this ratio 4>q is a few, so we adopt 0o = 2. 

6.2. Disk-Average Quantities 

Before considering the radial structure of the disk, let 
us consider the evolution of the galaxy as a whole be- 
tween z = 2 and 2 = 0. Our model does not allow the ro- 
tation curve or outer radius of the disk to evolve in time. 
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Fig. 1. — Time evolution from the beginning to the end of the 
fiducial simulation of the radially-integrated gas fraction, 2D Jeans 
mass at r = 8 kpc, and the radially-integrated star formation rate. 

However, over this redshift range, the circular velocity 
(assuming a constant spin parameter) will evolve by less 
than about 10% (e.g. iCacciato et aLll2011|) . Meanwhile, 
the position of the outer edge of the disk has a minimal 
effect on its evolution, so long as the outer edge of the 
star-forming disk is resolved. At larger radii than this, 
there is so little star formation that the gas is free to flow 
inwards at a constant rate and arrive at the edge of the 
star-forming disk unaltered by its passage through the 
HI disk. 

The primary changes in the disk are the steady de- 
cline in the accretion rate, and the steady formation of 
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stars. For the fiducial model, M ext (t) drops smoothly 
from about 13 M Q /j/r at z = 2 to 2 M Q /j/r at z = 0. 
This falloff is mirrored in the drop in total gas fraction, 
two-dimensional Jeans mass, and total star formation 
rate (figured]). The star formation rate in particular has 

almost the same numerical value as M ext (t), starting off 
slightly higher and converging to the accretion rate. This 
is a reflection of the fact that the formed stars can only 
come from gas that started in the simulation or accreted 
at a later time, and the initial gas reservoir is depleted 
in about 1 Gyr. 

Stars, once formed, remain in the disk, while the mass 
of gas in the disk falls with the cosmological accretion 
rate. This drives a steady decrease in the gas fraction 
from its initial value, down to 20%. Referring to the 
equilibrium solution for constant gas fraction (equations 
I5T)1 and H)D|) . and noting that f g has dropped by a factor 
of a few, while the accretion rate has dropped by a factor 
of about 6, we might expect a to decrease by maybe a 
factor of 2, while £ might decrease by more than a factor 
of 3. 

The two-dimensional Jeans mass ([Kim fe Ostrikerl 
l200l is defined by 



Gas 



Stars 



Mj 



G 2 E 



(64) 



Physically this represents the characteristic mass of a 
clump of gas which collapses under gravitational insta- 
bility to form a cluster of stars. Its steady decrease with 
time reflects the cooling of the disk, which allows smaller 
regions to collapse. This is the phenomenon that ex- 
plains why z ~ 2 galaxies contain giant clumps far larger 
than the biggest GMCs in present-day Milky Way-like 
galaxies. As a practical matter, this means that the typ- 
ical size of star clusters steadily decreases, so, coupled 
with the fact that a clump of gas can form stars with at 
most tens of percent efnency, clusters with M > 10 6 A/ Q 
are unable to form in today's quiescent spirals. In the 
fiducial model, Mj ~ 2 • 1O 7 M at r = 8 kpc. The de- 
crease in the upper envelope of cluster mass with time is 
consis tent with the arguments made bv lEscala fc Larsonl 
([2008]) . 

6.3. Radial Structure of the Disk 

We show the radial structure of our fiducial disk in fig- 
ures [21 E] and 01 We can understand the qualitative be- 
havior shown in these plots by considering the processes 
that drive the evolution. The two most important drivers 
are that Q = 1 almost everywhere at all times, and that 
stellar migration tends to self-regulate the stars such that 
Q* = Qum - recall that Qu m is a free parameter, below 
which stars are subject to transient spiral instabilities. If 
Q* > Qiim, stars will form and drive up £», decreasing 
Q*, while if Q* < Qu m , the stars will migrate inwards 
increasing er* and hence Q*. These two restrictions set 
Q g to a value somewhat less than Qu m , depending on 
the local ratio <7 9 /u*. These forces lead the simulations 
to form three qualitatively distinct regions: a stabilized 
stellar-dominated region, a star-forming region, and an 
HI disk. 

The radial extent of the star-forming region is more 
or less set by where the gas is molecular, i.e. fu 2 ~ 1- 
This in turn corresponds to where the gas column density 
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Fig. 2. — A direct comparison of the gas and stellar components 
as a function of radius at redshifts 2 (orange dotted), 1.5 (blue 
dot-dash), 1 (red dashed), and (black solid). The gas cools and 
depletes, while the stars accumulate and heat. The expanding 
stabilized region of the disk is evident in the dramatic decrease 
in gas transport velocity, large Q g , and a — > at. The outward 
movement of the region where stars form and migrate follows the 
peak in gas column density - Q* approaches Qu m = 2.5, the stellar 
metallicity gradient steepens, and the stellar scale height flattens. 

is larger than some metallicity-dcpendcnt critical value. 
For our fiducial initial conditions, the disk is molecular 
out to r w 15 kpc at z = 2. Within this radius, al- 
most the entire disk is vigorously star-forming. As time 
passes, a stellar-dominated central region begins to ap- 
pear. This occurs because, towards the center of the 
disk, the gas has short local dynamical times and hence 
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Fig. 3. — Radial profiles of quantities at redshift 2 (dotted), 1.5 
(dot-dashed), 1 (dashed), and (solid). The peak of fjj and hence 
the star formation rate move outwards as the simulation evolves, 
as the gas further in has been depleted and cannot be replenished. 

undergoes rapid star formation. In contrast, the inward 
mass flux of gas required to maintain Q w Qf is nearly 
independent of radius. Star formation depletes this gas 
as it moves inwards, so by the time it reaches the inner 
region of the disk, not only is there less gas than there 
would have been neglecting star formation, but it is be- 
ing consumed faster. In order to maintain a constant Q, 
given that ~ Qu m , the gas must maintain Q g close to 
constant. Star formation decreases the gas column den- 
sity, so to keep Q g roughly unchanged, the gas velocity 
dispersion must fall proportionally. Thus the gas velocity 



dispersion drops fastest in the center of the disk. 

By assuming a fixed gas temperature, we essentially set 
a floor on the value of a. When a hits this floor, which 
happens first at the inner edge of the computational do- 
main (see figure [2|) , the radiative loss rate C approaches 
zero. The gas no longer loses energy through shocks, and 
therefore ceases to move inwards. In this situation that 
region of the disk ceases to become gravitationally un- 
stable, and Q is allowed to rise. Without any means of 
mass transport, the gas simply depletes as it forms stars. 
As the gas column density drops off, the stars dominate 
the local stability of the disk. Since they are constrained 
to Q* ss Qum by our assumptions about stellar migra- 
tion, the overall value of Q/T of the disk in this region 
approaches Qu m as well. 

The third qualitatively distinct region of the disk may 
be thought of as the HI disk wherein fg 2 is low enough 
that stars form at a relatively slow rate, and gas flows 
in adhering even more closely to the equilibrium con- 
ditions of equations ([59]) and (|60|). which were derived 
by neglecting star formation in KB10, than in the star- 
forming region. In essence, the gas is allowed to flow in 
with a constant mass flux at each radius, since star for- 
mation is not depleting the gas significantly. Depending 
on the initial conditions of the simulation, the column 
density of stars may be low enough or the velocity dis- 
persion of the stars high enough that > Qu m for the 
duration of the simulation. In this situation the overall 
stability of the disk is almost exclusively determined by 
the stability of the gas, therefore the gas properties will 
correspond more closely to the equilibrium values with 
the gas fraction set to unity. 

Looking at the values for S and a near the solar circle 
(see figure [U, we see that they are too high relative to 
their observed values of approximately 13Mq/pc 2 and 8 
km/s respectively, though not by more than a factor of 
two. Moreover, the column density of gas near the center 
of the disk is lower than observed in the Milky Way. Both 
of these problems stem from the fact that when C — > 0, 
mass transport due to gravitational instability shuts off, 
whereas the real Milky Way has a number of mecha- 
nisms to transport gas into its central regions even when 
a at- The gas could be transported by a bar instabil- 
ity from larger radii, or the gas which we assume accretes 
at the edge of the disk could be accreting directly into 
the central region of the galaxy. Gas can also be recy- 
cled back to the ISM from stars. We assume this occurs 
instantaneously, so we neglect gas from stars which form 
farther out in the disk and migrate inwards. Nonethe- 
less, our model qualitatively reproduces the structure of 
z = disk galaxies: a central stellar-dominated bulge, an 
extended star-forming disk, and an outer Hl-dominatcd 
disk with very little star formation. 

6.4. Stellar Populations 

As the stars form in the fiducial simulation, one can 
treat them as adding together into a single population for 
the purposes of evaluating the torque equation, while at 
the same time evolving a number of passive populations, 
binned by age, alongside the single population. Only the 
active population affects the stability of the disk, while 
the passive populations simply serve as tracers of the 
stars formed during a particular epoch. This in turn is 
a reflection of the state of the gas at that time, with 
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Fig. 4. — Radial profiles of quantities at redshift 2 (dotted), 1.5 
(dot-dashed), 1 (dashed), and (solid). Within the star-forming 
region, the size of the Jeans mass decreases steadily, but increases 
at the center of the disk owing to the extremely low gas column 
densities. The two-component Q value transitions from unity in the 
gas (both H2 and HI) dominated regions to Q = Qu m T = 15/4 in 
the stellar-dominated component. 

the added effect of gradual stellar heating through radial 
migration. 

Stellar migration occurs locally as the result of star for- 
mation, since it is star formation which drives Q* below 
Qiim- It is therefore unsurprising that the stellar popu- 
lations seem to have very similar column density profiles 
(see figure [5]) to the star formation rate profile shown in 
figurc[3l The primary effect of migration is thus not mass 
transport inwards, so much as an increase in the veloc- 



ity dispersion. This can be quite significant - the oldest 
stars near the center of the disk reach nearly <r*.i = 100 
km/s, which is significantly larger than the gas velocity 
dispersion at any time in the simulation. 

The state of these populations near the solar neighbor- 
hood at z = is of particular interest, since these pop- 
ulations are well-observed and display well-known cor- 
relations. The velocity dispersions of stars in the so- 
lar neighborhood vary from about 17 km/s for 1 Gyr- 
old stars to ~ 10 Gyr - old stars with a* sa 37 km/s 
(|Nordstr5m et all 120041 : iHolmberg et~aT1 |2009T >. The 
th eoretical explanations for this correlation go back 
to iSpitzer fc Schwarzschildl (|1953l ) and generally center 
around the scattering of stars by molecular clouds and 
spiral structure, which gradually heats the disk. Other 
explanations h ave in cluded minor or maj or me rgers (e.g. 
Dierickx et"aDl2010t iBekki fc Tsuiimotoll2011t IQu et al.1 
201 lh and popping star clusters (|Assmann et al.l l20lTf iT 
All of these explanations are conceptually trying to do 
the same thing - form a thick disk from a thin disk. How- 
ever, a gravitationally unstable disk subject to star for- 
mation and a decreasing accretion rate will start with a 
high gas velocity dispersion that will decrease with time. 
This will also naturally generate an age- velocity disper- 
sion correlati on. This is the s cenar io presented in the 
simulations of lBournaud et all (12009 ), and in the chemo- 
dynamical models of iBurkert et al.l ( 19921 ). 

The age-velocity dispersion produced in our fiducial 
model may be explained as the combination of two phys- 
ical effects. First, the gas velocity dispersion decreases 
with time as the disk cools. This may be understood 
from the fact that if Q and Q* are self-regulated to con- 
stant values, then Q g must remain close to constant, and 
so if E decreases, so must a. As the gas cools, the stars 
it forms will be cooler than previous generations of stars, 
leading to an age-velocity dispersion correlation. The 
second effect is the heating of stars via transient spirals 
to maintain Q* = Qu m - Although this is a scattering 
process which heats stars over time, there is never a thin 
disk which gradually forms a thick disk. 

To better discern the importance of each of these ef- 
fects, we can compare the stars produced by the fidu- 
cial model to runs with certain effects artificially turned 
off. The high and low constant accretion rate mod- 
els shown in figure [6] have M e xt(t) — 12.3M©/yr and 
M ext (t) = 2.34Af Q /j/r respectively, corresponding to the 
accretion rates at the beginning and end of the fiducial 
simulation. For simulations where migration is turned 
off, we plot the properties of the stars at their epoch of 
formation, rather than their properties at z = 0. Thus 
the dynamical effects of migration as it affects the stabil- 
ity of the disk remain unchanged as compared with the 
fiducial simulation. Figure [S] shows explicitly that the 
age- velocity dispersion correlation is strongly affected by 
the accretion history and the presence of stellar heat- 
ing. All of the scenarios are able to generate some age- 
velocity dispersion correlation. Even the case with no 
stellar heating and a constant accretion rate produces 
one as the result of a fall in E, and hence a, as a result 
of star formation. 

7. DISCUSSION 

Starting from conservation laws and simple assump- 
tions about the gravitational stability of the disk, we 
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Fig. 5. — All stellar populations produced in the fiducial model 
at redshift zero, colored by their age with redder stars older. The 
ages are linearly spaced in time, so each population is about 1 Gyr 
of star formation. The dotted lines represent the initial population 
of stars, which has only evolved via stellar migration over the whole 
course of the simulation. Each newer population is less massive, 
dynamically colder, and has a steeper mctallicity gradient than its 
older analogues. 

have derived evolution equations for the radial profile 
of a two-component disk. Compared to semi-analytic 
models, this approach has the advantage that the vast 
variation in the state variables as a function of radius is 
resolved rather than averaged over the whole disk. This 
improvement comes with additional computational costs; 
however, these are not severe - even using the full Rafikov 
Q and multiple stellar populations, the code can evolve 
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Fig. 6. — Properties of stellar populations as a function of their 
age at a radius of 8 kpc at redshift zero. Note that the stars 
comprising the initial condition of the disk are not plotted here. 
Each line shows the result of a different model: the fiducial model 
(black), stellar migration off (red), high constant accretion (or- 
ange), low constant accretion (blue), stellar migration off and low 
constant accretion (purple). The models with constant accretion 
history are dashed. Every simulation produces an age- velocity dis- 
persion correlation via some combination of increasing cr» 4 of ex- 
isting stars or decreasing cr, which makes the younger stars dynam- 
ically cooler. 

a disk from z = 2 to z = on a single processo r in a few 
days, and using the IRomeo fc Wiegert] (|2011D approxi- 
mation to Q reduces the computation time to under an 
hour. 

This paper is primarily meant to introduce our 
methodology. However, the fiducial model demonstrates 
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a key point which is often overlooked in galaxy evolution 
and studies of the thick disk, namely that thick disks 
need not be formed from thin disks. An age-velocity 
dispersion correlation appears in our simulation, not be- 
cause of external perturbers, mergers, or gradual heating 
of a thin disk, but because a decreases with time and 
newly formed stars induce transie nt instabilities in the 
disk (see also iBurkert et al1l!992[ ). Both of these pro- 
cesses are strongly dependent on the cosmological situ- 
ation in which the disk finds itself, that is, its accretion 
history. Simulations of isolated thin disks which are grad- 
ually heated are therefore unrealistic, in the sense that 
they are missing the most important drivers of thick disk 
formation. The smooth increase in stellar velocity dis- 
persion with age produced in our simulations agrees qual- 
itatively with recent observations which demonstrate the 
lack of a distinctive bimodal ity between thick and thin 
disk stars ()Bovv et al.H2011aj) . 

This approach has several further applications which 
we intend to explore in future work. For Milky Way-like 
galaxies, even modern chemodynamical models with so- 
phisticated treatments of stellar migration and evolution 
rely on h ighly parameterized treatments of gas inflow in 
the d isk (jSchonrich fc Binneyl[200l iSpitoni fc Matteuccil 
1 2 1 If ) . If the gas evolves to keep the disk marginally grav- 
itationally unstable, its movement in the disk is not this 
simple - it depends on the evolution of the full non-linear 
set of equations we have derived here. By accounting for 
the diffusion of stars in radiu s as the result of scatterin g 
across corotation resonances (jSellwood fc Binnevl 12001) . 
our model could be extended to model the Milky Way 
in detail and compare directly with observations of the 
meta llicity gradient as a function of height above the 
disk dCheng et alJl2011|). the a ge- velocity dispersion cor- 
relation (jHolmberg et al.H2009D. the age-met allicity rela- 
tion or lack thereof (jEdvardsson et al.lll993f ).and the ra- 



dial an d vertical stellar density distributions (|Bovv et al.l 
1201 lbl) . 

Galaxy bimodality - the separation of galaxies into a 
blue cloud of star-forming galaxies and a red sequence of 
ellipticals - is often viewed as an evolutionary sequence. 
Blue cloud galaxies gradually accrete gas and smaller 
galaxies, which fuel star formation. Some of these galax- 
ies will undergo major mergers, leaving red and dead 
elliptical galaxies. These early-type galaxies can subse- 
quently undergo dry mergers, which extend the red se- 
quence to include e xtremely massive galax i es. B eyond 
this canonical view, ISanchez Almeida et al.l ()2011[ ) have 
noted the existence of a significant population of red spi- 
rals. 

By taking more realistic accretion histories from cos- 
mological simulations, we expect that a certain fraction 
of disks in the course of their lifetimes will experience a 
period of low accretion during which they will exhaust 
their gas supply and become redder, only to return to the 
blue cloud with the resumption of higher accretion rates. 

Given a set of realistic non-smooth but quiescent accre- 
tion histories, appropriate for a large fraction of sub-L* 
galaxies, we may therefore be able to reproduce aspects 
of phenomenology from the local universe out to z = 2 
as semi-analytic models do, but with the added benefit 
of a physical treatment of the disk dynamics. 
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APPENDIX 

NON-DIMENSIONAL EQUATIONS 

For the purposes of implementing the governing equations in a numerical code, it is useful to non-dimensionalize 
the equations. To do so is straightforward, and basically amounts to rescaling lengths to the radius of the disk, 
velocities to the circular velocity, and mass fluxes to the initial accretion rate of gas from the IGM. We can make 
the following substitutions, following KB10: r = xR, t = T[2TrR/v c j ) (R)], T = TM ext> ov^(R)R, dj = SjV<p(R), and 
= SjM ex t y o/(v<t,{R)R)- Here the subscript j may refer to gas or one of possibly many stellar populations. With 
these substitutions, the gas evolution equations (|4|) and ([7]) become 



dS_ 
df 
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dT 



(/3 2 + p + xj3'y - x(j3 + l)r" 



(fn. + M) 
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Primes denote partial derivatives with respect to x, and as with dimensional quantities, S and s with no subscript 
refer to properties of the gas. The dimensionless initial accretion rate is 
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GM e 



v^{R) 3 



(A3) 



The dimensionless thermal gas velocity dispersion is s t . 
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Employing the same procedure for the evolution equations of each stellar population's column density, we obtain 



= Jr -3=- + , (A4) 



dT J " \ dT J dT 



■w =87r n fH26ffKo T\l 1+ ^ (A5) 

AC . Mi 9 / ,/ OA 

-w =- 2 *»{ s '4 +s '-> + ir) (A6) 

where we have explicitly separated the effects of stellar migration and star formation. The dimensionless radial 
component of the bulk stellar velocity is y = v r */v^,(R). 
Similarly, the velocity dispersion evolution equations are 

f)a . 8s SF f)<j Mi 9 

~&T~~&T ~df~ ' ( ' 

-or * ,f «25ry-* ) -BT ' (A8) 

"OT = " 27r H^^ + S *J (A9) 

The change in velocity dispersion as a result of star formation is only an approximate relation, since it relies on a first 
order Taylor series expansion of the exact change in s*^, which in turn requires that S* t i ^> {dS*,i/dT) Mia dT. This 
condition cannot be satisfied when a completely new population of stars is formed as the simulation crosses into a new 
age bin, at which time S#_i = 0. Therefore we use the exact relation, 



l {S.,islJ old + f R (dSZf) S * 



where dSf J = dT(dS/dT) SF 
Finally we have the equations describing the transport of metals in the gas, 



dZ 2tt d\nZ , yu(l-f R ) dSj F 

dT (/3 + l)xSu dx T S dT [ ' 

and in a stellar population, 

!)7 Mig 

= -2nyS'^. (A12) 
The stellar metallicity change owing to the formation of new stars can be computed exactly as 

7 _ (SvjZ+^old + J R (dSf F )Z fA1Q\ 

>~>*,i,old ~r ttJ^ j 

These equations, given a torque r and a radial stellar velocity y, fully describe the evolution of the system. To obtain 
these two quantities, we imposed conditions on the evolution of Q and Q* (equations 1161 and I34[) . In dimensionless 
form these partial differential equations are 

yl | ( u 2 (l + /3) < | 51 | 1\_ max(Q Km -Q„0)M 

92 t" + 9i t' + g T = g F (A15) 
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where the coefficients of the dimensionlcss torque equation are 
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Both partial differential equations require an outer boundary condition, which essentially specifies the flux of each 
type of material at the edge of the disk. The mass flux of the gas is specified by some accretion history M ex t(i), 



r'(x = 1) 



Mextit) 



(l + /3(z = l)), 



(A20) 



while the flux of stars is set to zero via y(x = 1) = 0. The torque equation also requires an inner boundary condition, 
which we take to be t(x = xq) = 
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